home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Spearmans_Rho.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-09-13  |  14.3 KB  |  781 lines
  1. /* Spearmans Rho Test */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols>2 then do
  43.     'Message "Only two columns allowed!"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47.  
  48. /* Get cell reference for output range */
  49. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  50. if out_cell="" then do
  51.     'DEFPUBSCREEN("Workbench")'
  52.     exit
  53. end
  54. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  55.     'Message "Invalid cell reference"'
  56.     'DEFPUBSCREEN("Workbench")'
  57.     exit
  58. end
  59. /* Suppress Screen Redraw to Speed Things Up */
  60. 'Refresh 0'
  61.  
  62. /* Open a small output window on tcalc screen*/
  63. fo=0
  64. CR='0a'x
  65. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  66. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  67.      call writeln(6Info, DisplayMsg)
  68.     fo=1
  69. end
  70. else do
  71.     'Message "TCALC Screen not available for Progress messages"'
  72. end
  73. CALL DELAY(150)
  74.  
  75. /* Get cell references for top cell in each column */
  76. 'SelectCell' start_cell
  77. do col=start_col to end_col
  78.     'GetCursorPos'
  79.     top_cell.col=result
  80.     'Column 1'
  81. end
  82.  
  83. /* Get labels for later use on output */
  84. 'SelectCell' start_cell
  85. 'GetValue'
  86. testlabel=result
  87. testlabel=strip(testlabel)
  88. if datatype(testlabel,'n')=1 then do
  89.     labelflag=0
  90.     do x=1 to NCols
  91.     title.x="Column "||x
  92.     end
  93. end
  94. else do
  95.     labelflag=1
  96.     NRows=NRows-1
  97.     do x=1 to NCols
  98.     'GetValue'
  99.     str=result
  100.     title.x=translate(strip(str),"_"," ")
  101.     'Column 1'
  102.     end
  103. end
  104. if fo then call writech(6Info,"Progress...10 ")
  105.  
  106. /* Get data from cell range */
  107. col=start_col
  108. lav=0
  109. tot=0
  110. count.=0
  111. total.=0
  112. temp.=0
  113. do x=1 to NCols
  114.     'SelectCell' top_cell.col
  115.     if labelflag=1 then 'CursorDown 1'
  116.     do y=1 to NRows
  117.         'GetValue'
  118.         valtest=result
  119.         if datatype(valtest)='NUM' then do
  120.             'GetValue'
  121.             val=result
  122.             data.x.y=val
  123.             temp.x.y=val
  124.             count.x=1+count.x
  125.         end
  126.         'CursorDown 1'
  127.     end
  128.     col=col+1
  129.     val=0
  130. end
  131. if fo then call writech(6Info,"20 ")
  132.  
  133. /* Sort Temp Values */
  134. call Sort()
  135. if fo then call writech(6Info,"40 ")
  136.  
  137.  
  138. /* Calculate Ranks for Temp array */
  139.  
  140. Rank.=0
  141. NTies=0
  142. Do x=1 to NCols
  143.     N=count.x
  144.     Do y=1 to NRows
  145.     z=y+1
  146.     if temp.x.y=temp.x.z then do
  147.         cnt=1
  148.         beginrank=y
  149.         avrank=y
  150.         endrank=0
  151.         do until endrank~=0
  152.             y=y+1
  153.             avrank=avrank+y
  154.             z=y+1
  155.             cnt=cnt+1
  156.             if temp.x.y~=temp.x.z then endrank=y
  157.         end
  158.         if cnt>1 then NTies=NTies+1
  159.         avrank=(avrank/cnt)
  160.         do y=beginrank to endrank
  161.             Rank.x.y=trunc(avrank,2)
  162.         end
  163.     y=y-1
  164.     end
  165.     else do
  166.         Rank.x.y=y
  167.     end
  168.     end
  169. end
  170. if fo then call writech(6Info,"60 ")
  171.  
  172. /*Compare data to temp ranks and get data ranks */
  173. Rk.=0
  174. DO x=1 to NCols
  175.     Do y=1 to NRows
  176.         j=1
  177.         Do j=1 to NRows
  178.             if data.x.y=temp.x.j then do
  179.                 Rk.x.y=Rank.x.j
  180.                 leave j
  181.             end
  182.         end
  183.     end
  184. end
  185. if fo then call writech(6Info,"80 ")
  186.  
  187. /* Calculate Rho */
  188. TR=0
  189. SR=0
  190. SX=0
  191. SY=0
  192. SQ=0
  193. TT=0
  194. SD1=0
  195. SD2=0
  196. ZM=(NRows*((NRows+1)/2)**2)
  197. /*    DO y=1 to NRows
  198.         SX=(Rk.1.y)-((NRows+1)/2)
  199.         SY=(Rk.2.y)-((NRows+1)/2)
  200.         SQ=SQ+SX*SY
  201.     end
  202.     Rhoa=SQ/((NRows*((NRows**2)-1))/12)
  203. */
  204. if NTies~=0 then DO
  205.     in=2
  206.     de=1
  207.     Do x=1 to NCols
  208.         Do y=1 to NRows
  209.             total.x=(total.x)+Rk.x.y
  210.         end
  211.     end
  212.     /* Calculate Means */
  213.     mean.=0
  214.     do x=1 to NCols
  215.         mean.x=(total.x)/NRows
  216.     end
  217.     say mean.1
  218.     say mean.2
  219.     /* Calculate diffs */
  220.     diff.=0
  221.     sum.=0
  222.     dat=0
  223.     meenx=0
  224.     Do x=1 to NCols
  225.         meenx=mean.x
  226.         Do y= 1 to NRows
  227.             dat=Rk.x.y
  228.             diff.x.y=dat-meenx
  229.             sum.x=(dat-meenx)**2+(sum.x)
  230.         end
  231.     say sum.x
  232.     end
  233.     diffMult.=0
  234.     summulti=0
  235.     Do y=1 to NRows
  236.         diffMult.y=(diff.de.y)*(diff.in.y)
  237.         summulti=summulti+(diffMult.y)
  238.     end
  239.     rs=0
  240.     pr=0
  241.     rs=((summulti)**2)/((sum.in)*(sum.de))
  242.     Rhob=sqrt(rs)
  243. end
  244.     DO y=1 to NRows
  245.         TR=TR+((Rk.1.y)-Rk.2.y)**2
  246.     end
  247.     Rhoa=1-(6*TR)/(NRows*((NRows**2)-1))
  248.  
  249.  
  250. /* calculate Rho Probability */
  251. r=0
  252. ProbRho=0
  253. IF NTies=0 then do
  254.     ProbRho=1-PRHO(NRows,Rhoa)
  255.     if Rhoa~=1 then t=(Rhoa*sqrt(NRows-2))/(sqrt(1-Rhoa**2))
  256. end
  257. Else do
  258.     ProbRho=1-PRHO(NRows,Rhob)
  259.     if Rhob~=1 then t=(Rhob*sqrt(NRows-2))/(sqrt(1-Rhob**2))
  260. end
  261. df=NRows-2
  262. If Rhoa~=1 & Rhoab~=1 then do
  263.     /* Calculate T distribution function */
  264.  
  265.     P=PROBT(t,df)
  266.     if t>=0 then P=1-P
  267.     PT=P*2
  268.     /* Calculate T Critical */
  269.  
  270.     TCRIT1=TCRITICAL(.95,df)
  271.     TCRIT2=TCRITICAL(.99,df)
  272.     TCRIT3=TCRITICAL(.975,df)
  273.     if fo then call writech(6Info,"90 ")
  274.     TCRIT4=TCRITICAL(.995,df)
  275. end
  276. if fo then do
  277.     call writeln(6Info,"100")
  278.     call writeln(6Info,"Writing output to window...")
  279. end
  280.  
  281. /* Output */
  282. 'SelectCell' out_cell
  283. 'ColumnWidth 17'
  284. 'Put "Spearmans Rho - Rank Correlation Non-Parametric Test"'
  285. 'CursorDown 1'
  286. Do x=1 to NCols
  287.     title="Rank_"||title.x
  288.     'Alignment 2'
  289.     'Put' title
  290.     'CursorDown 1'
  291.     Do y=1 to NRows
  292.         'Put' Rk.x.y
  293.         'CursorDown 1'
  294.     end
  295. 'CursorUp' NRows+1
  296. 'Column 1'
  297. end
  298. 'SelectCell' out_cell
  299. 'CursorDown' NRows+3
  300. 'Put "Ties:"'
  301. 'Column 1'
  302. 'Put' NTies
  303. 'Column -1'
  304. 'CursorDown 1'
  305. 'Put "Rho (no Ties):"'
  306. 'Column 1'
  307. 'Put' format(Rhoa,,4)
  308. If NTies~=0 then do
  309.     'CursorDown 1'
  310.     'Put' format(Rhob,,4)
  311.     'Column -1'
  312.     'Put "Rho (corrected for ties):"'
  313.     'Column 1'
  314. end
  315. 'CursorDown 1'
  316. 'Put' format(ProbRho,,6)
  317. 'Column -1'
  318. 'Put "Prob.(Rho<=rho):"'
  319. If Rhoa~=1 & Rhoab~=1 then do
  320. 'CursorDown 1'
  321. 'GetCursorPos'
  322. new_cell=result
  323. 'Put "d.f.:"'
  324. 'CursorDown 1'
  325. 'Put "t:"'
  326. 'CursorDown 1'
  327. 'Put "P(T<=t) one-tail:"'
  328. 'CursorDown 1'
  329. 'Put "T-Critical (95%):"'
  330. 'CursorDown 1'
  331. 'Put "T-Critical (99%):"'
  332. 'CursorDown 1'
  333. 'Put "P(T<=t) two-tail:"'
  334. 'CursorDown 1'
  335. 'Put "T-Critical (95%):"'
  336. 'CursorDown 1'
  337. 'Put "T-Critical (99%):"'
  338. 'SelectCell' new_cell
  339. 'Column 1'
  340. 'Put' df
  341. 'CursorDown 1'
  342. 'Put' format(t,,4)
  343. 'CursorDown 1'
  344. 'Put' format(P,,6)
  345. 'CursorDown 1'
  346. 'Put' format(TCRIT1,,4)
  347. 'CursorDown 1'
  348. 'Put' format(TCRIT2,,4)
  349. 'CursorDown 1'
  350. 'Put' format(PT,,6)
  351. 'CursorDown 1'
  352. 'Put' format(TCRIT3,,4)
  353. 'CursorDown 1'
  354. 'Put' format(TCRIT4,,4)
  355. end
  356. 'Refresh 1'
  357. 'Refresh 2'
  358. /*'Message' "Finished"*/
  359. /*indicate the main script is finished*/
  360. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  361. result=writeln(6Info, DisplayMsg)
  362. if result~=0 then do
  363.     /*Wait 3 seconds*/
  364.     CALL DELAY(150)
  365.     /* close window*/
  366.     result=close(6Info)
  367. end
  368. 'DEFPUBSCREEN("Workbench")'
  369. exit
  370.  
  371. /* Procedures */
  372.  
  373. cellrow: procedure
  374. do
  375.     parse arg cell
  376.     do charpos=2 to length(cell)
  377.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  378.     end
  379.     return 0
  380. end
  381. Return
  382.  
  383. cellcol: procedure
  384. do
  385.     parse arg cell
  386.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  387.     cell=upper(cell)
  388.     len=length(cell)
  389.     val=0
  390. do charpos=1 to len
  391.     if datatype(substr(cell,charpos,1),n) then
  392.     do cell=reverse(substr(cell,1,charpos-1))
  393.     do x=1 to length(cell)
  394.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  395.     end
  396.     return val
  397.     end
  398.     end
  399.     return 0
  400. end
  401. Return
  402.  
  403. /* It is important to put the exposed array at the end of the next line */
  404. Sort: procedure expose NCols NRows temp.
  405. do x=1 to NCols
  406. L=(xtoy(2,int(log(NRows)/log(2))))-1
  407.     Do Until L<1
  408.     L=trunc(int(L/2))
  409.     Do J=1 to L
  410.             Do K=J+L To NRows By L
  411.             I=K
  412.             dumdat=temp.x.I
  413.             Do while I>L
  414.                 y=I-L
  415.                 If temp.x.y ~> dumdat then Leave
  416.                 temp.x.I=temp.x.y
  417.                 I=I-L
  418.             End
  419.             temp.x.I=dumdat
  420.             End
  421.         End
  422.     End
  423. End
  424. Return
  425.  
  426. syntax:
  427.      if arg(1)='FAIL' then do
  428.         'Message "Library is unavailable."'
  429.         'DEFPUBSCREEN("Workbench")'
  430.         exit
  431.         end
  432.     'DEFPUBSCREEN("Workbench")'
  433.     exit
  434.  
  435. Format:  procedure
  436.  
  437.      arg number, before, after
  438.      CallLine = SIGL
  439.      if ~datatype(CallLine, 'N') then CallLine = '??'
  440.  
  441.      /* Make sure we have a number as first (required) argument    */
  442.      if ~datatype(number, 'N') then do
  443.         if number = '' then
  444.            fc = 17     /* Wrong number of arguments           */
  445.         else
  446.            fc = 47     /* Arithmetic conversion error             */
  447.         signal FormatSyntaxError
  448.      end
  449.      num = number + 0
  450.      if before = '' & after = '' then
  451.         return num
  452.      else do
  453.         parse var num integer '.' fraction
  454.         if before = '' then before = length(integer)
  455.         if after = '' then after = length(fraction)
  456.         if ~datatype(before, N) | ~datatype(after, N) then
  457.            do fc = 18
  458.            signal FormatSyntaxError
  459.        end
  460.         if before < length(integer) then do
  461.            fc = 18
  462.            signal FormatSyntaxError
  463.         end
  464.         if after ~= length(fraction) then do
  465.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  466.         if integer<1&integer>-1 then integer=integer
  467.            else integer = integer + (fraction % 1)
  468.            fraction = substr(fraction, 3)
  469.         end
  470.         if fraction >= 0 then
  471.            return right(integer, before)'.'fraction
  472.         else
  473.            return right(integer, before)
  474.      end
  475.  
  476.  FormatSyntaxError:
  477.         if show('F', STDERR) then
  478.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  479.         else
  480.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  481.         'Message' mess
  482.         parse source Func .
  483.         if Func = 'FUNCTION' then do
  484.         'DEFPUBSCREEN("Workbench")'
  485.            exit "Err"
  486.         end
  487.         else do
  488.         'DEFPUBSCREEN("Workbench")'
  489.            exit 10
  490.         end
  491.  
  492. /*        Algorithm AS 89   Appl. Statist. (1975) Vol.24, No. 3, P377.
  493. c       
  494. c        To evaluate the probability of obtaining a value greater than or
  495. c        equal to is, where is=(n**3-n)*(1-r)/6, r=Spearman's rho and n
  496. c        must be greater than 1
  497. c
  498. c     Auxiliary function required: ALNORM = algorithm AS66
  499. */
  500.  
  501. PRHO: Procedure
  502.  
  503.     parse arg n, r
  504.     is=(n**3-n)*(1-r)/6
  505.     l.=0
  506.     zero=0
  507.     one=1
  508.     two=2
  509.     six=6
  510.     c1=0.2274
  511.     c2=0.2531
  512.     c3=0.1745
  513.     c4=0.0758
  514.     c5=0.1033
  515.     c6=0.3932
  516.     c7=0.0879
  517.     c8=0.0151
  518.     c9=0.0072
  519.     c10=0.0831
  520.     c11=0.0131 
  521.     c12=0.00046
  522. /*Test admissibility of arguments and initialize*/
  523.     prho = one
  524.     if (n <= 1) then return 0
  525.     if (is <= 0) then return 0
  526.     prho = zero
  527.     if (is > n * (n * n -1) / 3) then return 0
  528.     js = is
  529.     if (js ~= 2 * (js / 2)) then js = js + 1
  530.     if (n <= 6) then DO
  531. /* Exact evaluation of probability*/
  532.     nfac = 1
  533.     do i = 1 to n
  534.         nfac = nfac * i
  535.         l.i = i
  536.     end
  537.     prho = one / nfac
  538.     if (js ~= n * (n * n -1) / 3) then return
  539.     ifr = 0
  540.     do m = 1 to nfac
  541.         ise = 0
  542.         do i = 1 to n
  543.             ise = ise + (i - l.i) ** 2
  544.         end
  545.         if (js <= ise) then ifr = ifr + 1
  546.         n1 = n
  547.         DO UNTIL m=nfac
  548.             mt = l.1
  549.             nn = n1 - 1
  550.             do i = 1 to nn
  551.                 j=i+1
  552.                 l.i = l.j
  553.             end
  554.             l.n1 = mt
  555.             if (l.n1 ~= n1 | n1 = 2) then leave m
  556.             n1 = n1 - 1
  557.         end
  558.     end
  559.     prho = ifr / nfac
  560.     return prho
  561. end
  562. Else Do
  563.     b = one / n
  564.     x = (six * (js - one) * b / (one / (b * b) -one) -one) * sqrt(one / b - one)
  565.     y = x * x
  566.     u = x * b * (c1 + b * (c2 + c3 * b) + y * (-c4+ b * (c5 + c6 * b) - y * b * (c7 + c8 * b- y * (c9 - c10 * b + y * b * (c11 - c12 * y)))))
  567.     prho = u / exp(y / two) + znorm(x,1)
  568.     if (prho < zero) then prho = zero
  569.     if (prho > one) then prho = one
  570. end
  571. return prho
  572.  
  573. /*Algorithm AS66 Applied Statistics (1973) vol22 no.3*/
  574.  
  575. /*Evaluates the tail area of the standardised normal curve
  576.   from x to infinity if upper is .true. or
  577.   from minus infinity to x if upper is .false.*/
  578.  
  579. znorm: procedure
  580.     parse arg x, upper
  581.     zero=0
  582.     one=1
  583.     half=.5
  584.     ltone=7.0    
  585.     utzero=18.66
  586.     con=1.28
  587.     p=0.398942280444
  588.     q=0.39990348504
  589.     r=0.398942280385
  590.     a1=5.75885480458
  591.     a2=2.62433121679
  592.     a3=5.92885724438  
  593.     b1=-29.8213557807
  594.     b2=48.6959930692
  595.     c1=-3.8052e-8
  596.     c2=3.98064794e-4
  597.     c3=-0.151679116635
  598.     c4=4.8385912808
  599.     c5=0.742380924027
  600.     c6=3.99019417011  
  601.     d1=1.00000615302
  602.     d2=1.98615381364
  603.     d3=5.29330324926 
  604.     d4=-15.1508972451
  605.     d5=30.789933034
  606.       up=upper
  607.       z=x
  608.     if (z<zero) then do
  609.         up=0
  610.               z=-z
  611.     end
  612.        if (z>ltone | up=0 & z>utzero) then Do
  613.               alnorm=zero
  614.               if (up=0) then alnorm=one-alnorm
  615.         return alnorm
  616.     end
  617.     y=half*z*z
  618.     if (z<=con) then do
  619.         alnorm=half-z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3))))
  620.         if (up=0) then alnorm=one-alnorm
  621.         return alnorm
  622.     end
  623.     alnorm=r*exp(-y)/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6))))))
  624.     if (up=0) then alnorm=one-alnorm
  625. return alnorm
  626.  
  627. PROBT: PROCEDURE
  628.  
  629.     PARSE ARG TA,K1
  630.     A=.36338023
  631.     W=ATAN(TA/SQRT(K1))
  632.     S=SIN(W)
  633.     C=COS(W)
  634.     L=K1-2*INT(K1/2)
  635.     IF L=0 THEN DO
  636.         T1=S
  637.         IF K1=2 THEN DO
  638.         PO=.5*(1+T1)
  639.         RETURN PO
  640.         END
  641.         T2=S
  642.         J1=-1
  643.         J2=0
  644.         K2=(K1-2)/2
  645.     END
  646.     ELSE DO
  647.         T1=W
  648.         IF K1=1 THEN DO
  649.             T1=T1*(1-A*L)
  650.             PO=.5*(1+T1)
  651.             RETURN PO
  652.         END
  653.         T2=S*C
  654.         T1=T1+T2
  655.         IF K1=3 THEN DO
  656.             T1=T1*(1-A*L)
  657.             PO=.5*(1+T1)
  658.             RETURN PO
  659.         END
  660.         J1=0
  661.         J2=1
  662.         K2=(K1-3)/2
  663.     END
  664.     DO I=1 TO K2
  665.         J1=J1+2
  666.         J2=J2+2
  667.         T2=T2*C*C*J1/J2
  668.         T1=T1+T2
  669.     END
  670.     T1=T1*(1-A*L)
  671.     PO=.5*(1+T1)
  672. RETURN PO
  673.  
  674. TCRITICAL: PROCEDURE
  675.  
  676.     PARSE ARG PO,K1
  677.     AO=.0001
  678.     E=.005
  679.     E2=E+E
  680.     A=2*PO-1
  681.     IF K1=1 THEN DO
  682.         TO=TAN(1.57079633*A)
  683.         RETURN TO
  684.     END
  685.     IF K1=2 THEN DO
  686.         SN=SIGN(2/(1-A*A))
  687.         IF SN=-1 THEN TO=-1*(A*SQRT(ABS(2/(1-A*A))))
  688.         ELSE TO=A*SQRT(2/(1-A*A))
  689.         RETURN TO
  690.     END
  691.     P1=PO
  692.     Z=NORMALPP(PO)
  693.     W=Z*(1+2/(1+8*K1))
  694.     T3=K1*(EXP(W*W/K1)-1)
  695.     SELECT
  696.         WHEN Z<0 THEN TT=-1
  697.         WHEN Z=0 THEN TT=0
  698.     OTHERWISE TT=1
  699.     END
  700.     SN=SIGN(T3)
  701.     IF SN=-1 THEN T3=-1*(TT*SQRT(ABS(T3)))
  702.     ELSE T3=TT*SQRT(T3)
  703.     /* LABELA:
  704.     TO=T3
  705.     PO=PROBT(TO,K1)
  706.     F1=PO-P1
  707.     TO=T3+E
  708.     PO=PROBT(TO,K1)
  709.     F2=PO
  710.     TO=T3-E
  711.     PO=PROBT(TO,K1)
  712.     F2=F2-PO
  713.     F2=F2/E2
  714.     T4=T3-F1/F2
  715.     IF ABS(T4-T3)>AO THEN DO
  716.         T3=T4
  717.         SIGNAL 'LABELA'
  718.     END
  719.     */
  720.     T4=0
  721.     DO FOREVER
  722.         TO=T3
  723.         PO=PROBT(TO,K1)
  724.         F1=PO-P1
  725.         TO=T3+E
  726.         PO=PROBT(TO,K1)
  727.         F2=PO
  728.         TO=T3-E
  729.         PO=PROBT(TO,K1)
  730.         F2=F2-PO
  731.         F2=F2/E2
  732.         T4=T3-F1/F2
  733.         IF ABS(T4-T3)>AO THEN T3=T4
  734.         ELSE LEAVE
  735.     END
  736.     TO=T4
  737.     
  738. RETURN TO
  739.  
  740. LOGGAMMA: PROCEDURE
  741.  
  742.     ARG XA
  743.     C1=76.18009173
  744.     C2=-86.50532033
  745.     C3=24.01409822
  746.     C4=-1.231739516
  747.     C5=.001208580
  748.     C6=-.000005364
  749.     C7=2.506628275
  750.     X1=XA-1
  751.     W=X1+5.5
  752.     W=(X1+.5)*LN(W)-W
  753.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  754.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  755.     L=W+LN(C7*S)
  756. RETURN L
  757.  
  758. NORMALPP: PROCEDURE
  759.  
  760.     ARG P0
  761.     A1=2.515517
  762.     A2=.802853
  763.     A3=.010328
  764.     A4=1.432788
  765.     A5=.189269
  766.     A6=.001308
  767.     Q=.5-ABS(P0-.5)
  768.     SN=SIGN(-2*LN(Q))
  769.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  770.     ELSE W=SQRT(-2*LN(Q))
  771.     W1=A1+W*(A2+A3*W)
  772.     W2=1+W*(A4+W*(A5+A6*W))
  773.     ZZ=W-W1/W2
  774.     SELECT
  775.         WHEN (P0-.5)<0 THEN TT=-1
  776.         WHEN (P0-.5)=0 THEN TT=0
  777.         otherwise TT=1
  778.     END
  779.     ZZ=ZZ*TT
  780. RETURN ZZ
  781.